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In the past year we have made significant progress in improving our fundamental under- 
standing of the physics of this problem, as detailed below. Furthermore, having brought our 
code to a state of fairly robust functionality, we devoted significant effort to optimizing it for 
running long cases. We optimized the code for vectorization to the extent that it now runs 
eight times faster than before (a typical case used to take a substantial fraction of a Cray 2 
hour to run to convergence). 

Physical improvements to the model 

In starting to model very dense particle layers, when the particle mass density can exceed 
100 times the local gas mass density, we realized that in such regions, the viscosity arising 
from interparticle collisions may become comparable to, or even exceed, the turbulent gas 
viscosity, and began to explore realistic implementations of particle viscosity v v . One simple 
parametrization of v v is the particle viscosity routinely used in planetary rings which are not 
overly optically thick ( e.g . Goldreich and Tremaine 1978, Wisdom and Tremaine 1988): 
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where v p is the average particle random velocity, Cl is the orbital frequency, and / is the vertically 
integrated particle area filling factor or optical depth. The trick is to estimate v p . A simple 
ring-type assumption such as v p = CIS, where 6 is the particle layer thickness, is inappropriate 
in this case since global particle motions on the scales 6 and Cl may be driven by turbulent 
eddies without nearby particles having much relative velocity at all if they are sufficiently well 
coupled to the local gas velocity. 

We have come up with a simple model for v p based on the particle stopping time t p , which 
determines the Schmidt number 5c.The model begins with the scaling relationships 

v p ~< V ><1 >~< u >< l > 2 ~< v > 2 / < u >, (2) 

where <«>,</ >, and < lj > axe the characteristic relative velocity, length scale, and collision 
frequency of the particles respectively. There are two regimes of interest. If the collision time 
and length scales are not limited by the global scales of the system (layer thickness 6 and orbital 
frequency Cl), then u ~ nirr 2 v Te i ~ v re j//*, where l* is the mean free path between collisions. 
From the work of Volk et al (1980), we identify v rt i as a fraction \Jtpftg < 1 of the global 
average particle random velocity v p when t p <t g , where t p is the particle stopping time and t g 
is the eddy turnover time. Naturally, v re ; can never exceed v p , and we account for this in the 
limit t p » t g . We also note from our Schmidt number model (cf. also Volk et al 1980) that 
( v fl/ t, p) 2 = Sc = (1 + t p /t g ), or v p ~ v g Sc~b. Consequently, v p = l*v Tt i = l*v p T{Sc), where 

T(Sc) = ' /SC ~ 1 if t p < t g , and :F(Sc) = if t p > t g . (3) 
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Note that this implies the particle viscosity goes to zero for very small t p regardless of the 
particle density. This is because the particles are all trapped to the same local gas velocity and 
have no relative random velocity at all. 

In another limit, if the scales of the system limit the collision frequency or mean free path 
(as in, for instance, optically thin planetary ring systems), v p ~ / < v >< l > or with the same 
assumptions as above, v p = fv t /Sc where v t is the turbulent gas viscosity. 

We have implemented viscous terms in our numerical code using these parametrizations of 
v p . For the cases we have been studying (30 - 100 cm radius particles, 1 and 10 AU, minimum 
mass solar nebula) the particle optical depth / is on the order of 0.1 and the Schmidt number 
is on the order of unity; consequently the particle viscosity is about 10% of the gas turbulent 
viscosity. 

Numerical results 

Using the newly vectorized code, we ran several models which included particle viscosity 
terms both at 1 AU and at 10 AU. The code is well behaved in both limits. The results differ 
from previous runs in that mean radial and azimuthal velocities in the particle layer are now 
more slowly varying with vertical distance from the midplane due to the increased coupling by 
particle viscosity. 

Reynolds averaging of fluid equations 

One aspect of our model that we wanted to put on firmer ground is our mixed use of Favre 
(mass) averaging for the momentum equations with Reynolds (time) averaging for the particle 
conservation equation. It is the Reynolds averaging that results in the diffusion term we use 
to model particle layer diffusion, whereas similar terms are suppressed in the Favre averaged 
equations. Feeling that this was rigorously inconsistent and perhaps even quantitatively im- 
portant, we have devoted considerable effort to rewriting the momentum equations from the 
Reynolds- average standpoint. At this time we have obtained the new correlation terms but not 
as yet coded them up. They are all tractable and can be modeled in very much the same way as 
the gas eddy viscosity is always derived as a model of the Reynolds stresses, and heat transport 
and particle diffusion by convection are modeled with the gradient diffusion hypothesis (using 
the Prandtl and Schmidt numbers respectively). The terms are small, and we expect them to 
change our numerical results in minor but potentially very interesting ways. 

Modeling of turbulence and viscosity 

In our prior work, we have explored two independent parametrizations of the nebula turbu- 
lence, of differing complexity. In Champney and Cuzzi (1990), we pointed out the poorer than 
desirable agreement between the eddy viscosity as determined from the two-equation ( k — e) 
model and that obtained with our current Prandtl model, which is characterized by only one 
parameter (the critical Reynolds number). Since the two-equation model is partly ad hoc and 
contains at least five constants, we have chosen so far to use the simpler Prandtl model. How- 
ever, as we pointed out last year, not even the Prandtl model is without its uncertainties. The 
critical Reynolds number Re* (Champney and Cuzzi 1990, equation 49) depends on the nature 
of the flow regime. Heretofore we have used a value of 500 for Re*, but now believe that the 
true value of Re* is about 100. Use of this number brings the two-equation and Prandtl models 
into agreement. 

However, the Prandtl technique cannot model the damping of turbulence by the particle 
phase (Sproull 1961, Elghobashi and Abou-Arab 1973, Pourahmadi and Humphrey 1983); this 
may be very important not only in the shear layer, but also in earlier phases of the nebula 
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when particle settling and accretion occurs in the presence of widespread convective turbulence. 
Consequently we are delving more deeply into self-consistent turbulence models. 

The two-equation models currently in use (e.a. Rodi 1984) postulate one equation for the 
generation, transport, and damping of the turbulent kinetic energy k, and a similar equation 
for the energy dissipation rate e. We have verified that the Jb-equation (including its particle 
damping terms) is derivable in a straightforward way from the basic fluid equations, while as 
far as we can determine, the e-equation is a relatively ad hoc creation designed to improve fits 
to data. In fact, prior to the current widespread use of k-c models, use of only the ^-equation 
was standard (Rodi 1984). We find that the dissipation rates e calculated with our two-equation 
model are approximated by k/t g , where the eddy turnover time t g is simply the inverse of the 
orbital frequency. In the coming year, we plan to replace the e-equation entirely by the simple 
scaling c = k/t g in the ^-equation, simplifying the method to a one-equation model. This will 
expedite the study of turbulence in the shear layer, including particle damping. 
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